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Abstract 


The response of the living tissue to the effects of strong heating or 
cooling can cause the blood flow rate to vary by an order of magnitude. A 
mathematical model for the freezing of living tissue is formulated which 
takes into account the nonlocal temperature dependence of the blood 
flow rate when the temperature distribution in the tissue is substantially 
nonuniform, as in cryosurgery. 


INTRODUCTION 

Mathematical analysis and prediction of temperature distribution in living 
tissue during the process of freezing has been used in the study and optimization 
of cryosurgical procedures. In the last years a number of different approaches 
to describing the heat transfer process in living tissue have been proposed (for 
a review see|Q-|15 ). Within the framework of these approaches the obtained 
bioheat equation is of the 


dT 

CPUT = V(fcVT) - CbPbJfiT -Ta) + S , 
at 


( 1 ) 


where T is the tissue temperature, c and p denote the specific heat and den¬ 
sity of the tissue, Cf, and pb are the specific heat and density of blood, J is 
the blood flow rate per unit tissue evolume, (i.e. the volume of blood flowing 
through a unit volume of tissue per unit time), k is the thermal conductivity of 
the tissue Ta is the systemic arterial blood temperature and S is the 

rate of metabolic heat generation. Cofactor /, ranging from 0 to 1 is due to 
heat exchange between arterial and venous blood flowing through the nearest 
vessels ||5|-[p^. 
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Propagation of the freezing front P is conventionally described in terms of 
the free boundary problem of the Stefan-type ^ : 


-(fcV„r) |r+ +(fcV„T) |r_ , 

(2) 

T\r^=T\r_=Tf , 

(3) 


where L is the latent heat of fusion, r+ and P- denote the boundaries of the 
freezing front on the living and frozen sides of the tissue, respectively, and Tf 
is the freezing temperature (fig.l). 

A more accurate description of the heat transfer process in living tissue 
is obtained if one takes into account the fact that living tissue form an active, 
highly heterogeneous medium. In particular, the phase transition in living tissue 
does not occur at a single temperature, but over a temperature range. The 
influence of this effect has been investigated in [^j. Also, when the size of the 
frozen region of the tissue is small in comparison with the characteristic length 
of the blood vessels that directly control the heat exchange between the cellular 
tissue and blood, the heterogeneity of living tissue has a substantial effect on 
the heat transfer process. Therefore, equation (1) which models the blood flow 
rate in terms of a continuous field J(r) has to be modified. 

The response of the living tissue to the effects of strong cooling or strong 
heating can cause the blood flow rate to vary by an order of magnitude Q. 
In general, if the temperature distribution in the tissue is substantially nonuni¬ 
form, as for example in cryosurgery, then the temperature dependence of the 
blood flow rate is nonlocal and the blood flow rate at a given point depends 
on certain characteristics of the temperature distribution rate, other than the 
tissue temperature at the point only. Indeed, variations in the blood flow rate at 
a given point are caused mainly by the overall response of the blood collection 
vessels, involving arteries of different length and radius. Hence, equation (1) 
should be modified to include a relationship specifying the tissue temperature 
dependence on the blood flow rate § . 

In this paper we shall formulate a phenomenological model for the thermal 
response in living tissue during the freezing process that will include the effects of 
the aforementioned factors. It should be pointed out that the main equations to 
be formulated can be obtained by rigorous analysis of the microscopic relations 
governing the heat transfer process in living tissue and the response of the 
vascular network, which we shall consider in a subsequent paper. 

TWO BOUNDARY MODEL 

Let us assume that all characteristic scales of the temperature distribution 
over the living tissue are greater than the mean length ly of the vessels that 
directly control the heat exchange between the cellular tissue and blood. In 
this case equation (I) is actually the result of averaging over the scale ly the 
microscopic equation describing heat exchange between blood in different vessels 
and the cellular tissue. Therefore, the quantity J appearing in equation (1) is 
actually the blood flow rate averaged on the scale ly , rather than the true blood 
flow rate Jt. 
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Figure 1: Two boundary model for tissue freezing processes: a - the vessel 
resistance as a function of the blood temperature, b - characteristic region of 
living tissue. 


Nonuniformity in the true blood flow rate Jt can be explicitly associated 
with the distortion of the vascular network and, therefore, can be characterized 
by scales much less than ly. On scales substantially larger than ly the averaged 
blood flow rate and the true blood flow rate coincide, leading to the identity J = 
Jt- On scales less than ly, the averaged blood flow rate obeys the conservation 
law for the total blood flow. 

A simplest relation that allows for the features mentioned above can be 
expressed in the form|5[|-|p^ 

J-V(/2VJ) = Jt . (4) 

The form of the second term on the left-hand side of (4) follows from the re¬ 
quirement that its integral over the whole tissue domain under consideration 
should be equal to zero. The mean length ly of the vessels controlling the heat 
transfer between blood and the cellular tissue in a given domain of radius ly is 
determined by the mean blood flow in this domain, i.e., by the averaged blood 
flow rate J. 

To find the relationship between ly and J, let us consider the heat transfer 
in a given domain of diameter 1. We assume that such a domain, on the average, 
contains one artery (arteriole) and one vein (venule) of length I, which is typical 
is the case of living tissues . The characteristic rate of heat generation (or 
dissipation) in a domain due to convective blood flow in these vessels can be 
estimated as CbPb{T* — Ta)Tra^v, where T* and v are the mean temperature and 
velocity of blood in the vessels and a is their radius. Similarly, the characteristic 
rate of heat generation due to heat conduction in living tissue is given approx¬ 
imately by cp{T — Ta)V ■ k/{pcP), where 1/ ~ is the volume of the domain, 
and Pcp/k is the characteristic time of temperature variation on the scale I, 
caused by heat conduction. On the scale Ic the two quantities should be of the 
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same order, i.e. 


cp{T* 



CbPb(T - Ta)Tra‘^v 


and the blood temperature in the corresponding vessels must be in thermody¬ 
namic equilibrium with the cellular tissue, i.e., T* T. Taking into account 
that J ^ Tra^vjVc and cp ~ CbPb we obtain 


I 


2 

c 


A i 

pc J 


(5) 


In particular, for typical values of kjpc ~ 2 • 10“^ cm^/s and J ^ 3 • 10“^ s“^ 
it follows from (5) that Zc ~ 1 cm. 

Given the true blood flow rate Jt, (equation (1),(4)) and the relation (5) 
determine the temperature distribution in living tissue. In order to complete 
the model we must describe the response of the vascular network to the tissue 
temperature variation. Since the heat exchange between cellular tissue and 
blood in the vessels of length larger than Ic is negligibly small , the convective 
heat flow through such vessels obeys the conservation law. Therefore, the heat 
flow CbPb(T* — To) li (li = 7r(a^v)i) in a vein i of length I > Ic is approximately 
equal to the total heat dissipation rate due to blood flow in the domain 
which is directly supplied with blood through this vein and the corresponding 
artery. In other words, for such a vein i we have 


CbPb {T* - Ta) Ii= J dr CbPb J (T - To) , 

Qi 

where the averaged blood flow rate J can be replaced by the true rate Jt, since 
I ^ ly. Upon substitution, this equation is valid formally for veins whose length 
I <ly. Indeed, in this case it reduces to the identity 

/i = y Jtdr (6) 

Qi 

and blood in such a veins is in the thermodynamic equilibrium with cellular 
tissue, i.e. T* ~ T. Thus, for any vein i we may write 

{T* - Ta) h = J dr Jt (T - Ta) . (7) 

Qi 

A valid response of the vascular network to temperature variation requires 
that the resistance of each artery i should depend on the mean tissue tem¬ 
perature in the domain rather than on particular details of temperature 
distribution over the domain Qi. Taking into account equation (7), it is natural 
to assume that an artery i responds to the temperature T* of the blood in a vein 
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i, because the value of T* can be treated as a mean tissue temperature in the 
corresponding domain Qi. In other words, at least in the quasistationary case, 
we may suppose that the resistance Ri to blood flow in a vein i is an explicit 
function of (T* — To). 

The blood flow distribution {7^} over the vascular network is directly deter¬ 
mined by the vascular network architectonics and total pressure drop across it, 
provided all the vessel resistances are have given. Thus, when the distribution 
of the temperature dissipation rate Jt {T — To) in the living tissue is a given 
held, the explicit dependence of the blood temperature T* in the vein i on the 
blood how li (through the dependence of the vessel resistances on the blood 
temperature, see (7)) completely determines the blood how pattern {Ii\. In 
other words, the temperature dissipation rate can be treated as the information 
held for the thermoregulation of living tissue. 

Hence, in the general case the evolution equation for the true blood how rate 
may be represented as 

r> T 

T^ + Jt = Jo + F{JtiT-Ta)}, (8) 

where the transient term allows for a possible time delay r in the vascular 
network response to the temperature variations, Jq is a uniform blood how rate 
when T = Ta and F{Jt{T—Ta)'\ is certain functional operator of the information 
held Jt{T — Ta) that specihes thermoregulation in the living tissue. The simplest 
expression for the response operator F, which models thermoregulation correctly 
may be written in the form 

F{Jt{T - To)] = 1 / dr'G(r,r')Jt(r')[r(r') - T,] , (9) 

Qo 

where Qo is the total living tissue domain under consideration, G(r,r') is the 
kernel of the operator F, r and r' are vectors, and A is the length of temper¬ 
ature survival of the living tissue, in particular A « Ta — Tf. Typically, the 
temperature response of the vessels within a tumor is strongly depressed. This 
allows us to set 

G(r,r')=0, (10) 

where the point r' belongs to the tumor domain Qt,i.e. r' € Qt- When the 
tissue temperature attains a certain low value T^r > Tf, the vessels exhaust 
their ability to respond to temperature variation and the blood flow rate is no 
longer dependent on the tissue temperature. In mathematical terms this effect 
can be taken into account if in the response operator F{Jt{T — Ta)} the true 
tissue temperature T is replaced by a seeming tissue temperature Tg, which 
coincides with the true tissue temperature T when T > T„r and is equal to 
Tyr for Tf < T < Tyy. Thus, for the given linear response operator (9) the 
temperature response of the vascular network may be described by the equation 
with saturation by the following equation 

= Jo + ^ y Jr'G(r,r')Jt(r')[T,(r') - TJ , (11) 

Qo 


5 



where 



T if T > Tyr 
Tyr if Tf <T < Tyr 


( 12 ) 


In the frozen region Qf where T <Tf 


Jt = 0 . 


(13) 


The system of equations (1),(4),(11) with boundary conditions (2),(3) and rela¬ 
tions (12),(13) form the desired complete phenomenological model for the freez¬ 
ing process of living tissue. It should be pointed out that, within the framework 
of the proposed model, the living tissue freezing process may be represented 
by the propagation of two boundaries, which divide the tissue into three re¬ 
gions. The first one is the frozen region, which is separated from the extremely 
cooled tissue region (the second region) by the interface T where T = Tf. In 
the extremely cooled region the temperature self-regulation is depressed and an 
effective boundary of seeming temperature separates it from the third region, 
where the temperature varies signihcantly in the presence of high blood flow 
rate. 

To complete the description of our model for the freezing processes of living 
tissue we need to specify the kernel G'(r, r'). The kernel should accounts for 
two different effect. The first of these is increase in the blood flow rate at the 
point r, associated with additional blood flow coming from systematic circu¬ 
lation through the host artery of the vascular network. The second effect is 
the redistribution of the available blood flow in the vascular network. However, 
when the resistance to blood flow in the vascular network is determined by a 
group of vessels that very greatly in length, then the thermoregulation processes 
will give rise to the formation of a low resistance path on the vascular network 
that connects the host artery and the small arteries located in the region of 
considerable temperature increase. In this case the blood flow redistribution ef¬ 
fect is ignorable in comparison with first effect mentioned above, and the kernel 
(^(rjr') may be represented approximately by the d-function 


G(r,r') « 5(r - r')0+(r') . 


(14) 


where 0+ = 0 if r S Qt and 0+ = 1 for r ^ Qt- 

By definition, the thermal self-regulation process, where the kernel G(r, r') 
is represented by the expression (14), is called ideal thermoregulation. 

For ideal thermoregulation, equation (11) takes the following form for the 


unfrozen region 



(15) 


when T < Tyy and r ^ Qt, for r e Qt 


Jt = Jo ■ 


(16) 
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In the extremely cooled region, where Tf <T < Tyr and for r ^ Qt 


Jt — Jo 


Tf-Tg 

Tf - Tyr 


(17) 


Thus, for ideal thermoregulation the description of the freezing process for 
living tissue is reduced to the system of equation (1),(4),(15), the boundary 
conditions (2),(3) and the expressions (14),(16),(17). It should be noted that in 
this model the frozen region and the region where the blood flow rate strongly 
depends on temperature are separated by a layer (the extremely cooled region) 
where the blood flow rate has the maximum and is independent of the tissue 
temperature. 


References 

[1] R.I.Andrushkiw. Mathematical modelling of freezing front propagation in 
biological tissue. Math.Comp.Modelling V.B. NIO, 1-9, (1990). 


[2] G.Comini and S.Del. Gindice. Thermal aspects of cryosurgery. 
J.Heat.Transfer. 534-549, (1976). 

[3] T.E.Cooper and G.J.Trezek, Mathematical predictions of cryogenic lesion, 
in Cryogenics in Surgery, H. von Leden and W.G.Cohan (eds.), Medical 
Examination Publishing Co., New York, 128-149, (1971). 

[4] A.A.Gage, K.Guest, M.Montes, J.A.Caruana, D.A.Whalen, Jr., Effect of 
varying freezing and thawing rates in experimental cryosurgery, Cryobiol¬ 
ogy 22, 175-182, (1985). 

[5] I.Lubashevskii, V.Gafiychuk, A.Cajan. Bioheat Transfer. Book (to be 
published). 

[6] I.Lubashevskii, V.Gafiychuk. Mathematical modelling of heat and mass 
transfer in living tissue. In.: Proceedings of the IMACS symposium on 
Mathematical modelling. February 1994.Vienna. 

[7] I.Lubashevsky, A.Priezzev, V.Gafiychuk, M.Gadjan. Local thermal coagu¬ 
lation due to laser-tissue interaction as irreversible phase transition. Journal 
of Biomedical optics. 2(1), pp. 95-105. 1997. 

[8] I.Lubashevsky, A.Priezzev, V.Gafiychuk, M.Gadjan. Dynamic free bound¬ 
ary model for laser thermal tissue coagulation in.: Laser tissue interaction 
and tissue optics 11. G.Delacretar, Ed. Proc.SPIE. 2923 (1996) pp. 48-57. 


7 



[9] I.Lubashevsky, A.Priezzev, V.Gafiychuk.Free boundary model for local 
thermal coagulation. Growth of spherical and cylindrical necrosis domain. 
In: Laser tissue interaction VII. S.L.Jacques, Ed. Proc.SPIE. 2681 (1996) 
pp. 81-91. 

[10] AndrushkiWjR.L, Gafiychuk,V., Lubashevskii,!. Two boundary model for 
freezing processes in living tissue.in Proceedings of the 14th IMAGS World 
Gongress ,March 1994.Atlanta.pp.123-126. 

[11] G.I.Mchedlishvili. Microcirculation of blood general principles of control 
and disturbances. (Nayka. Leningrad 1989) (In Russian). 


[12] A.Shitzer and R.G.Eberhart, in Heat transfer in medicine and biology. 
Anal, and Appl. edited by A.Shitzer and R.G.Eberhart (Plenum, New- 
York, 1985) v.l, p.l37. 


[13] A.Solomon, On the use of mathematical models in cryosurgery. Computer 
Science Div., Oak Ridge National Laboratory, TN, (1980). 


[14] S.Weinbaum and L.M.Jiji. A new simplified bioheat equation for the 
effect of blood flow on local average tissue temperature. ASMI Journal of 
Biomechanical Engineering. v.l07, 131-139, (1985). 


[15] G.H.Wissler. Comments on the new bioheat equation proposed by 
Weinbaum and Jiji. ASMI Journal of Biomechanical Engineering. v.l09, 
226-233, (1987). 



